import os
os.environ["PROJ_LIB"] = r'C:\Users\ckpark\anaconda3\Library\share'

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import numpy.ma as ma
import scipy.stats as ss

def init_plotting():
    plt.rcParams['font.family'] = 'Times New Roman'
init_plotting()

def moving(inte, n):
	result = np.convolve(inte, np.ones((n,))/n, mode='same')
	return result

ycomp = np.arange(0,140.1,0.1)


######## climatology ##########
csyear = 1992
ceyear = 2021
dayy = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]

read = np.loadtxt('../DATA/BIN/47156', dtype='str')
yr = np.array(read[:,0], dtype='int')
mo = np.array(read[:,1], dtype='int')
dy = np.array(read[:,2], dtype='int')
pre = np.array(read[:,3], dtype='float')

'''
for t in range(len(yr)):
    if yr[t] == csyear and mo[t] == 1 and dy[t] == 1:
        stp = t
    if yr[t] == ceyear and mo[t] == 12 and dy[t] == 31:
        etp = t
        
dpre = pre[stp:etp]
cpre = np.zeros((365))
for t in range(365):
    cpre[t] = np.mean(dpre[t::365])
    
climd = []; cc = 0
for m in range(12):
    for d in range(dayy[m]):
        mo = m + 1
        dy = d + 1
        if mo == 12 and dy == 31:
            climd.append(cpre[cc])
        cc += 1
climd = []; cc = 0
for m in range(12):
    for d in range(dayy[m]):
        mo = m + 1
        dy = d + 1
        if mo >= 1:
            climd.append(cpre[cc])
        cc += 1
'''        
################################

read = np.loadtxt('edi_exp_input', dtype='str')
datas = np.array(read[:,:], dtype='float')
    

cdat = []; a= []
for t in range(len(datas[:,0])):
    for k in range(len(datas[0,:])):
        if datas[t,k] > 0.0:
            cdat.append(ycomp[k])
            a.append(datas[t,k])
            break

ddat = []; b = []
for t in range(len(datas[:,0])):
    for k in range(len(datas[0,:])):
        if datas[t,k] > 1.0:
            ddat.append(ycomp[k])
            b.append(datas[t,k])
            break

acdat = []; bcdat = []
for t in range(len(cdat)):
    if (t+1) == 31:
        print(t+1,cdat[t]*(t),ddat[t]*(t))
    if (t+1) == (31+28):
        print(t+1,cdat[t]*(t),ddat[t]*(t))
    if (t+1) == (31+28+31):
        print(t+1,cdat[t]*(t),ddat[t]*(t))        
        
    acdat.append(cdat[t]*(t+1))
    bcdat.append(ddat[t]*(t+1))

rcdat = []; rddat = []
for t in range(len(cdat)):
    rcdat.append(cdat[t]*10)
    rddat.append(ddat[t]*10)

fig = plt.figure(figsize=(10,10))
outfig = '../99.FIG/EDI_EXI_ideal.png'
for xx in range(1):
    ax1 = fig.add_subplot(1,1,1)
    #---------------------------
    xax = np.arange(0,len(datas[:,0]))   
    jet = plt.cm.get_cmap('Spectral')
    #jet = plt.cm.get_cmap('coolwarm_r')
    
    cmin = -2.5
    cmax = -1.0 * cmin  
    levels = np.linspace(cmin,cmax,9)
    
    #ax1.set_ylabel(ylabels[xx],rotation=1,fontsize=20,labelpad=90)
    #ax1.yaxis.set_label_coords(-0.11,0.2)
    
    ax1.set_yticks(np.arange(0,len(ycomp),100))
    ax1.set_yticklabels(ycomp[::100])
    ax1.set_ylim(0,1210)
    ax1.set_ylabel('Precipiation (mm)',fontsize=26)
    
    xticks = []
    xlabs = []
    xticks.append(0); xlabs.append(1)
    cc = 0
    for k in range(len(xax)-1):
        if ((k + 1) % 10) == 0:
            xticks.append(k)
            xlabs.append(k+1)
   
    ax1.set_xticks(xticks)
    ax1.set_xticklabels(xlabs)
    ax1.set_xlim(-0.5,89)
    ax1.set_xlabel('Acculumated date',fontsize=26)
    
    ax1.tick_params(labelsize=25)
    
    #datas = np.ma.masked_greater(datas, 2.5)
    ax2 = ax1.pcolormesh(xax-0.5,np.arange(0,len(datas[0,:]))-0.5,np.transpose(datas),cmap=jet,vmin=cmin,vmax=cmax)
    
    ax1.plot(xax,rcdat,color='black',linestyle='solid',linewidth=3)
    ax1.plot(xax,rddat,color='black',linestyle='dashed',linewidth=3)


    cbaxes = fig.add_axes([0.91, 0.14, 0.02, 0.7]) # left, bottom, width, height
    cb = fig.colorbar(ax2, cax=cbaxes,ticks=levels,orientation='vertical',extend='both')
    #cb = fig.colorbar(ax23,ticks=levels,fraction=0.040,pad=0.06)
    cb.ax.tick_params(labelsize=25)
    plt.subplots_adjust(hspace=0.0,wspace=0.2)
    plt.show()      
    fig.savefig(outfig,bbox_inches='tight')    


fig = plt.figure(figsize=(10,5))
outfig = '../99.FIG/EDI_EXI_accum_prec.png'
for xx in range(1):
    ax1 = fig.add_subplot(1,1,1)
    #---------------------------
    xax = np.arange(0,len(acdat))   
    
    #ax1.set_yticks(np.arange(0,5))
    #ax1.set_yticklabels(np.arange(0,5))
    ax1.set_ylim(75,400)
    ax1.set_ylabel('Accumulated precipiation (mm)',fontsize=26)

    xticks = []
    xlabs = []
    xticks.append(0); xlabs.append(1)
    cc = 0
    for k in range(len(xax)-1):
        if ((k+1) % 10) == 0:
            xticks.append(k+1)
            xlabs.append(k+1)
    
    ax1.set_xticks(xticks)
    ax1.set_xticklabels(xlabs)
    ax1.set_xlim(0,90)
    ax1.set_xlabel('Accumulated date',fontsize=26)
    
    ax1.tick_params(labelsize=25)
    
    ax1.plot(xax+0.5,bcdat,color='blue',linewidth=3,alpha=1.0,label='1.0 scEDI threshold')
    ax1.plot(xax+0.5,acdat,color='green',linewidth=3,alpha=1.0,label='0.0 scEDI threshold')
    
    #ax1.bar(xax,bcdat,color='gray',alpha=0.5,label='1.0 EDI threshold')
    #ax1.bar(xax,acdat,color='black',alpha=0.5,label='0.0 EDI threshold')
    
    plt.legend(fontsize=20)

    plt.subplots_adjust(hspace=0.0,wspace=0.2)
    plt.show()      
    fig.savefig(outfig,bbox_inches='tight')    

print(acdat[0],bcdat[0])